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A two-dimensional superconductor (SC) on surfaces of topological insulators (TIs) is a mix- 
ture of s-wave and helical p-wave components when induced by s-wave interactions, since spin 
and momentum are correlated. On the basis of the Abrikosov-Gor'kov theory, we reveal that 
unconventional SCs on the surfaces of TIs are stable against time-reversal symmetric (TRS) 
impurities within a region of small impurity concentration. Moreover, we analyze the stability 
of the SC on the surfaces of TIs against impurities beyond the perturbation theory by solving 
the real-space Bogoliubov-de Gennes equation for an effective tight-binding model of a TI. We 
find that the SC is stable against strong TRS impurities. The behaviors of bound states around 
an impurity suggest that the SC on the surfaces of TIs is not a topological SC. 

KEYWORDS: topological insulator, helical Dirac electron, unconventional superconductivity, impurity scat- 
tering, time reversal symmetry, spin orbit interaction, Abrikosov-Gor'kov theory, Bogoliubov- 
de Gennes equation, impurity induced state 



1. Introduction 

Three-dimensional (3D) topological insulators (TIs) 
have two-dimensional (2D) surface states (SSs) topologi- 
cally protected by the time reversal symmetry (TRS). 1,2 
The existence of SSs is characterized by Z2 topological 
invariants, which are determined by the band structure of 
the bulk. 3-6 In most cases, strong spin orbit interactions 
(SOIs) play important roles in constructing topologically 
nontrivial band structures and induce TIs. By SOIs, he- 
lical spin structures in momentum space are observed, 
i.e., helical Dirac electrons are generated. Helical Dirac 
electrons in TIs have been verified through observations 
of energy dispersions of Dirac SSs by angle-resolved pho- 
toemission spectroscopy (ARPES). 7 ~ 10 

Superconductivity on the surfaces of TIs attracts at- 
tention as one type of 2D unconventional superconduc- 
tivity. An unconventional superconductor (SC) is a mix- 
ture of s-wave and helical p-wave components, since in 
helical Dirac electron systems, spin and momentum are 
correlated. 11 Such a SC is possibly induced by a prox- 
imity effect from an s-wave SC to a TI 12 ' 13 or, in other 
words, by an s-wave attractive interaction, 11 when the 
Fermi energy is away from the Dirac point. This SC for- 
mally resembles a spinless chiral p-wave SC that breaks 
TRS 14 in the representation where the Dirac electron dis- 
persion is diagonalizcd. The difference is that the SC on 
the surfaces of TIs docs not break TRS. From that simi- 
larity of the two SCs, an unconventional SC on the SSs of 
TI is proposed for application to quantum computations 
using Majorana bound states caused by the proximity 
effect between a superconductor and the surface states 
of TI. 12 Introducing superconductivity into the surfaces 
of TIs has been a challenge in experimental research. For 
example, toward this goal, a SC has been realized in Cu- 
doped Bi2Sc3, 15: 16 although the existence of a surface SC 
is not confirmed yet. Recently, it has been reported that 
superconductivity is introduced into Bi 2 Se 3 thin films by 
the superconductivity proximity effect. 17 

In realizing such a SC, the stability of the SC is an 



important problem. In particular, impurity effects are 
relevant to the stability since surfaces on TI frequently 
contain disorders such as defects or impurity potentials. 
Moreover, on surfaces of TIs, there are facets or steps 
with disordered boundaries. For example, by scanning 
tunneling spectroscopy (STS) studies of the surfaces of 
Bi2Te3, one of TIs, an abrupt change in the local density 
of states (LDOS) was observed near a step structure. 18 
Moreover, studies of the stability against impurities 
allow us to clarify the fundamental physics of a SC on 
the surfaces of TIs since the stability of the SC depends 
on the symmetries of the order parameter and on impu- 
rities. For example, conventional s-wave SCs are robust 
to TRS impurities because pair breaking does not ex- 
ist if the impurity concentration is small. 19 In contrast, 
unconventional anisotropic SCs are fragile against TRS 
impurities because of anisotropic quasiparticle scatter- 
ing that induces pair breaking. 20-23 Such pair breaking 
is also induced even in s-wave SCs when impurities break 
TRS. 24 For the present SC, the order parameter is a mix- 
ture of s-wave and helical p-wave components. We study 
quasiparticle scattering in an unconventional SC in or- 
der to reveal its similarity to or difference from those of 
other SCs. 

The organization of this paper is as follows: In §2, in 
order to study the fundamental stability of the SC on 
the surfaces of TIs, we analyze the impurity concentra- 
tion dependence of the mean-field order parameter in the 
small concentration range in an idealistic helical Dirac 
electron model using a perturbation theory referred to 
as the Abrikosov-Gor'kov (AG) theory. 24 We find that 
such a SC, as well as conventional s-wave SCs, is robust 
in that the mean-field critical temperature T c and the 
mean-field order parameter An do not linearly decrease 
with the TRS impurity concentration. 25 

In §3, we investigate impurity effects nonperturba- 
tively by solving the real-space Bogoliubov-de Gennes 
(BdG) equation for the tight-binding model of Bi 2 Se 3 , 
which is an effective model of TI. 26,27 We study the scat- 
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tering strength and impurity concentration dependences 
of the SC. Moreover, we show that induced bound states 
around impurities are not Andreev bound states, which 
implies that the present SC is not a topological SC. We 
summarize our study in §4 with a discussion. 

2. Pert ur bat ive Approach 

In order to study impurity effects on the surfaces of 
TIs using the AG theory, we introduce a 2D helical Dirac 
electron dispersion as an effective model of SSs with an 
s-wave attractive interaction and on-site TRS scatter- 
ing following a previous letter. 25 We introduce the s- 
wave attractive interaction because such an interaction 
is the most well-known origin of superconductivity, for 
example, an electron-phonon interaction. Moreover, on- 
site TRS scattering is introduced because it is one of the 
simplest impurity scatterings and is useful to study fun- 
damental impurity effects. Thus, our Hamiltonian con- 
sists of three parts, i.e., a 2D helical Dirac electron dis- 
persion, Ho, an s-wave attractive interaction term, 'Hint, 
and an on-site TRS impurity scattering term, H lmp : 



H — Ho + H n 



Hu 



(2.1) 



In this section, we assume that the 2D helical Dirac 
electron dispersion is represented as 



Ho = ^2 c t (fc)[wF(o' K fc :r 

k 



(Tyky) 



»I]c(k), (2.2) 



with the Fermi velocity vp > 0. We set the Fermi en- 
ergy fj, > in order to consider the branch of Dirac 
electrons above the Dirac point (called the "+" branch 
hereafter) and neglect mixture of the two branches. Here, 
the Pauli matrix er describes the electron spin and c(k) = 

(Cfet Ckif '. 

We introduce the representation in the helicity ba- 
sis (helicity representation) to diagonalize Ho- By 
using the unitary transformation d/ kT = (cj^ + 
rexp(i#fc)cL )/y/2 (r = + or — ), Ho is diagonalized as 



k 



eft(k)(v F \k\T g -tt)d(k), 



(2.3) 



where t z is the z component of the Pauli matrix describ- 
ing branches of Dirac electrons, 9k = axg{k x + ik y ), and 
d(k) = (dk+ dk~) T ■ The index r (r = ±) represents 
branches of Dirac electrons. Then, we define the energy 
£fc as the energy of the "+" branch measured from the 
Fermi energy as £ fc = vp\k\ — \x. We neglect the "-" branch 
and write d? as d^ below. 

Note that the argument in this section is also appli- 
cable to the 2D helical Dirac electron dispersions with 
different spin-momentum relations. For example, in the 
SSs of Bi2Se3, one of TIs, a x k y — a y k x is substituted for 
Cxk x + &yk y in eq. (2.2). 26 ' 27 In such a case, we have to 
redefine 9k as 9k = axg{k x + \k y ) + |. 

Wc assume that the s-wave attractive interaction H mt 
is written as 

Hint = ^ E ^nt(fe,fe';S,s') C -fe s c L' C fe's' C -fe'^- 4 ) 



k.k' ,s.s f 



we assume that 

V int (k,k';iA) = V int (k,k'-A,i) 

v int (k,k'-AA) = v int (k,k'u,i) 



-g e [-w c ,' 

(6,,£fc' £ [-w c ,' 



-g' g hw c , 

[-w c , 



with a cutoff w c [i and g, g' > 0. Note that, under the 
condition uj c <^ ^, the interaction affects only electrons 
on the "+" branch. Then the interaction term in the 
helicity representation is 

jl_ 

4S 



Hi 



-j-Y^^'-^di k d{d k ^ 



(2.7) 



k,k' 



Here, we introduce an on-site TRS impurity scattering 
term as 

^nP = iEE ei(fc " fc) ^ ct ( fc ) C ( fc '), (2-8) 
i=l k.k' 

where N\ is the number of the impurities in the system 
and Ri (i = 1, • • ■ , Ni) is the impurity location. Then the 
impurity scattering term in the helicity representation is 



Ni 



Hirap g 



7? E E e-^P(e k 9 k+q )di +q d k , (2.9) 

i—1 k,q 

where P(9) = exp(i#/2) cos(0/2) is a phase factor specific 
to Dirac electron systems. Here, P(n) = means that the 
backscattering is forbidden. Moreover, this phase factor 
contributes to the it Berry phase, which leads to an an- 
tilocalization effect of single Dirac cone systems, 28 i.e., 
the electric conductivity in the system has a positive 
quantum correction. 

We introduce a mean- field approximation in eq. (2.7) 
and construct a BCS mean-field Hamiltonian. By the 
BCS-type decoupling, a momentum-dependent pair po- 
tential A(fe) is derived as 

■9 „-i0 fa 



A(fe) 



25 



k' 



e^'(d k ,d_k>) = Ae- lWfc , (2.10) 



and our interaction term is approximated as 
1 



Hi 



2E[ A ( fc ) rft - 



A*(fe)dfcd_ fc . (2.11) 



The pair potential in eq. (2.10) resembles a spinless chiral 
p-wave pair potential. This resemblance is related to the 
emergence of Majorana bound states around integer vor- 
tices. 12 ' 14 The difference between the two SCs is that the 
SC of helical Dirac electrons is time-reversal-symmetric, 
but the spinless p-wave SC is not. 

In the representation of the original electrons opera- 
tor c, this SC is composed of a mixture of s-wave (sin- 
glet) and p-wave (triplet) symmetries, and the interac- 
tion term is represented as 



Hi 



where S is the size of the 2D system. In this equation, 



E E K 

k S1.S2 



h.c. 



,(2.12) 
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where 



A(fc) = A s (fe) + A t (fc) 



A t (fc) 



A 
~2 



1 

-1 







(2.13) 



Here, the s-wave pairing A s and the p-wave pairing A t 
are introduced. The mixture of two components origi- 
nates from the broken inversion symmetry in the surface. 
The d vector corresponding to the triplet pairing compo- 
nent At has a momentum dependence d(k) oc (k x , k y , 0), 
i.e., At has a helical p-wave symmetry. Here, the helical 
p-wave SC is TRS p-wave SC, in which the chirality of 
the pair potential is different for each spin component. 29 
By using the mean-field approximation in eq. (2.11), 
the Hamiltonian is approximated as 

1 



H 



MF "I- V inlp I * 



2 ^ 

fc 



&(k)H MF (k)V(k) 



^EE ei(fc " fc)R,vI/t ( fc )^( fc ' fc ')*(fc')- 



25 



l — l fc,/c' 

(2.14) 

We introduce a Nambu representation 'l'(fc) = 

(^k,e ^k,h) T = (dk d\k) T ■ The ma t r i x elements in eq. 
(2.14) are defined as 

H MF (k) ■■ 



& Ae" 16 " 
A*e w <° 



Vi(k,k') = 







(2.15) 



We define an imaginary time Green's function as 
G(t) = - (T T *(r)* T (0)), where the time evolution 
of the operators is obtained from ^^\t) = 

e r ^ij/(t) e -T"'H_ Then the Green's function in frequency 
space is G(iu> n ) = dr e 1 "" T G(r). From the equation 
of motion of the Green's function, the Gor'kov equation 



(iu n - H)G(icj n ) = 1 



(2.16) 



is derived, where w n = (2n + 1)tvT is the fcrmionic Mat- 
subara frequency and T is the temperature. 

When we assume that the scattering term Vi mp is per- 
turbation using a perturbation series expansion with re- 
spect to u/S, the Green's function G is represented as 



G(iu} n ) = G (iuJn) 



(G (i0Jn)Vimp) n G (iuJn) , 



(2.17) 



where the nonperturbative Green's function Go(iw n ) 



(iu> n — Hmf) 1 is introduced. 

According to the AG theory, we perform the impurity 
average operation so that the system recovers its trans- 
lational symmetry. By this operation, diagonal terms 
of momentum remain in the Hamiltonian and we only 
have to consider diagonal terms in the Green's func- 
tion. By the impurity average operation, the quantity 
IS j=i e 1 ( q ' Ri ~ q '' R ^ is replaced with ?^i<5q ; q', where 
rij = Ni/S is the impurity concentration. We perform 
the impurity average of the right-hand side of eq. (2.17) 
for each term. The terms that correspond to n — 1 and 
n = 2 in eq. (2.17) arc represented as diagrams in Figs. 
1(a) and 1(b), respectively. They are the lowest-order 
terms about ni. However, the term shown in Fig. 1(a) is 
negligible for the estimation of T c because it just causes 
a constant self-energy shift and does not contribute to 
relaxation processes due to pair breaking. 



a 



(b) 




Fig. 1. Diagrams for self-energy terms of the lowest order in eq. 
(2.17). (a) and (b) correspond to n = 1 and n = 2, respectively. 
Wavy lines represent impurity scatterings. 



Therefore, we consider the diagram in Fig. 1(b) and 
its higher-order series as the self-energy term for the es- 
timation of T c . The contribution of this term is 0{rii). We 
adopt this approximation following Abrikosov-Gor'kov 24 
and Ambcgaokar-Griffin. 30 Then the Green's function is 
calculated as 

G(k,iuj n ) = 



G {k,iuj n ) 
1 



-±(k, 



iui n - HuF(k,ioj n ) 
where the self-energy term is 



(2.18) 



£(fe,iu;„) 



Vx (fc, k')G(k', iu n )Vi (k', k) 



fliU m—^ \Jr 
fc' n 



-o k )Y 



-A* n e ie " iu) n -£ k , 
and the renormalized mean-field Hamiltonian is 

Cfe A n e- ie »" 



H M F(k,iuj r , 



(2.19) 



(2.20) 



In eqs. (2.18) -(2.20), Cj n and A„ are the renormalized 
frequency and pair potential, respectively 

By using eqs. (2.18) -(2.20), ui n and A„ are self- 
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consistently calculated as 



1 



A n = A 



2r 
1 



I A, 



A, 



2r 2 



|A„| 2 



(2.21) 



(2.22) 



where t\ and r 2 are two types of relaxation times. In 
the case of the present SC, the two relaxation times are 
calculated as 

r2n d9 



1 



1 



2n 2t 2 



irriiU 2 No 



2tt 



cos 2 (6»/2) 



(2.23) 



where N = j, J2k <K£*) i s tri e density of states (DOS) 
at the Fermi energy. Therefore n and t 2 are the same 
when considering O(nj). 

In order to estimate the order parameter or critical 
temperature, we need a self-consistent equation. A self- 
consistent equation is obtained using eq. (2.10): 



A = - 



2/35 



5^e**G (!ifc (fc > iu B ). 



(2.24) 



The self-consistent equation for the order parameter 
at zero temperature, Aq, is 



Ao(rii) 



A 



g f Uc duj ^ 

2PS J_ Uc 2^ \ £ 2 +£ 2 +A 2 

and A = A 



(2.25) 



where ui = lu + ^ ; — 

On the other hand, the self-consistent equation for the 
mean-field critical temperature T c is 



gT c (rii) y-vy-v (2T 2 |6Dn|) 1 



2,5" 



(2.26) 



Because the long-range SC order does not develop in 2D 
systems at a finite temperature, T c calculated from the 
mean-field theory provides a criterion of the Berezinskii- 
Kosterlitz-Thoulcss (BKT) transition for the develop- 
ment of a quasi-long-range order. 31,32 

By using the above equations, the values are obtained 

as 

7r 2e< ( 2 

T c (nO = T c (0) - 



4r s (n;) 



, T c (0) = w c exp 



.9^0 



the results of the present SC in comparison with other 2D 
SCs. For each SC, (2ti)~ 1 and (2t 2 ) _1 are calculated. We 
call SCs with 7-7 1 = 0(nf) stable SCs, while SCs with 
(r,)" 1 = 0{n\) are fragile. In the table, "magnetic scat- 
tering" indicates that the impurity Hamiltonian takes a 
form as 



u 



Ni 

u x 

~ > c 



i=l 



(2.30) 



which represents the scattering by magnetic impurities 
polarized along the z-direction. The forms of matrix ele- 
ments for the BCS mean-field Hamiltonian HMp(k) and 
the scattering Hamiltonian Vi(k,k') are also shown in 
the table. 



(2.27) 



Ao(ni) - A (0) 



4r s (ni) 



, A (0) = 2w c cxp — 



2 

(2.28) 



where the relaxation time t s is defined by 

(r s )- 1 = (2r 1 )- 1 -(2r 2 )- 1 . (2.29) 

Because (r s ) _1 = 0(nf) is satisfied, T c and Ao do not 
decrease linearly in the present SC. 

We compare this result with the results for other full- 
gap SCs reported in the literature 23,24 . Table I shows 
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Table I. Stability of SCs: (a) 5-wave SC with TRS impurities, (b) 5-wave SC with magnetic impurities, (c) Chiral 
p-wave SC with TRS impurities, (d) SC on TIs with TRS impurities, (e) SC on TIs with magnetic impurities. 







Vi(k,k>) 


stability 


(a) 


U(k) a\ 
\A* -£ k ) 


io-ij 


oGdDIc 


l D J 


\a* -ay 


Olj 


II dgllc 


(r) 


\A*e»> -a ) 




II dgllC 


(d) 


( i{k) Ae~ Wk \ 


\ o -p(0_ fc -0_ fc o/ 


stable 


(e) 


/ £(fc) Ae~*M 


/ - #fc + 7T) \ 

\ -P(0_ fe - 0_ fe , + it)) 


fragile 
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To summarize this section, unconventional SCs in- 
duced by the s-wave attractive interaction on the sur- 
faces of TIs are robust to TRS impurities. This result is 
achieved by calculating the dependences of T c and A on 
the TRS impurity concentration, where T c provides a cri- 
terion of the BKT transition for a quasi-long-range order 
in 2D systems. The unconventional SC on the surface of 
TI is robust because of the cancellation of two phase fac- 
tors, one from the pairing potential and the other arising 
when a Dirac electron is scattered by a TRS impurity. 
In contrast, unconventional SCs reported in the litera- 
ture, such as the <i-wave and chiral p-wave SCs, 20-23 are 
sensitively suppressed through scattering by a tiny con- 
centration of impurities because of the phase factor of 
the pairing potential. 

We treated impurities as perturbations in this sec- 
tion. The perturbation theory is valid if u 2 n\N is much 
smaller than the pair potential A though a rough esti- 
mate. This estimation is derived from the reduction in 
the pair potential due to magnetic scattering. 

3. Nonperturbative Approach 

3.1 Model and method 

In order to study the impurity effects on the SC on the 
surfaces of TIs by real-space BdG calculation, we use a 
tight-binding model of Bi2Se3 in slab geometry, which 
is an effective model of 3D TI. 26 ' 27 We obtain a single 
Dirac cone on the surface of the slab. 

In order to study impurity effects, wc consider a Hamil- 
tonian composed of three terms in the same way as in 
52: 



7i — !~Lq + Tin 



(3.1) 



Here, Ho is the effective tight-binding Hamiltonian of 
Bi2Se3 , which is first introduced by Zhang et al. 26 On 
the other hand, H- ln t is the interaction term and H lmp is 
the impurity potential. 

The Hamiltonian Ho has a structure of a 4 x 4 ma- 
trix, because of the presence of two orbitals and spin 
indices. The two orbitals are the antibonding and bond- 
ing orbitals constructed from the p z -orbitals of Bi and 
Se atoms; they have different parities from each other. 
We refer to these orbitals as U E" and orbitals, re- 
spectively. Bi2Sc3 has a rhombohcdral symmetry, but for 
simplicity we adopt a model reduced to the D4h symme- 
try. 26 In studying fundamental properties of SSs near the 
Dirac point, this approximation is valid since anisotropy 
owing to the rhombohcdral symmetry is weak near the 
Dirac point. 

A Bloch representation of Ho 26 is 

H (k) 

(e k + M k A z k z A n k- \ 
A z k z e k -M k A\\k- 
A\\k+ e k + M k -A z k z 
\ A\\k + -A z k z e k -M k J 



(3.2) 



where 



<k 
M k 



D z k z , 



M-B\\k\-B z kl (B||,B,<0), 



k± — k x i iky, fc|| — {k x , ky), 

where A\\ and A z represent the strength of SOIs, while 
M is the energy difference between the two orbitals. The 
band curvatures of the two subbands arc different, which 
arises from nonzero B^ and B z . Here, we note that, due 
to the gauge transformation defined in eq. (18) of ref. 27, 
Ho{k) is not invariant under the operation of a standard 
choice for the time reversal operator, / ® io~ y ■ JC, where 
/ is the identity matrix acting on subband indices, o~ y is 
the Pauli matrix acting on the spin indices, and K, is the 
complex conjugate operator. 

In order to construct a Wannier representation of Ho, 
we perform a substitution such that 

ki — > sin ki , fcj — > 2 — 2 cos ki (i = x, y, z), 

and the Fourier transformation 



N x N v N~ 



c(fc) = 



c f (fc) = 



£££e^c(r), (3.3) 



N v N~ 



EEE^^)- ( 3.4) 



H a = ^c t (fc)i? (fc)c(fc), c(k) = (c fctE Cfc tH c fcJ , E C fc4 .H 



For simplicity, we choose the model parameters in eq. 
(3.2) as 

A x = 1 cV, A z = 0.5 eV, B x = 1.5 cV, B z = 0.3 cV, 

D x = D z =0,M = 0.5 eV, (3.5) 

where the particle-hole symmetry about subbands E and 
H is imposed. We assume that the lattice parameter is 
normalized as 1. Then the curvatures of two subbands 
are equivalent except signs. For readers, we cite the pa- 
rameters given by ab initio calculations 26 below: 

Ax = 0.5 eV, A z = 0.2 eV, B x = 0.6 cV, B z = 0.1 eV, 

D x = 0.1 cV, D z = 0.05 eV, M = 0.3 eV. (3.6) 

We obtain a single Dirac cone at the T point in a sur- 
face by considering a slab geometry when the param- 
eters support nontrivial topological indices. In this 
paper, we impose an open boundary condition in the z- 
direction and periodic boundary conditions in the x- and 
y-directions. Figure 2 shows energy dispersions for % - 
In the dispersion, gapless SSs exist inside the bulk band 
gap. 

Then, we introduce the s-wave attractive interaction 
term Hint as 

iijiCTjCr' ,t,t' 

where g T a T a ,{i,j) = g^ltiJ) holds. 

As the third term of our model Hamiltonian, we intro- 
duce the s-wave TRS impurity scattering Hamiltonian 
, T 'Himp- The impurity scattering Hamiltonian Hi m p is rep- 
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Fig. 2. Dispersion when the set of Hamiltonian parameters is 
A x = 1 eV, A z = 0.5 eV, B x = 1.5 eV, B z = 0.3 eV, D x = D z = 
and M = 0.5 eV. An open-boundary condition is imposed in 



resented as 



Ni 



"^imp — C l(i) C n(i) 



JJ 1 



imp 



C, 



(3.8) 



where N\ is the number of impurities, n(i) is the location 
of the i-th impurity, u is the strength of the impurity 
potential, and Himp is the matrix representation of the 
Hamiltonian. 

By solving the BdG equation for the introduced model 
Hamiltonian, we analyze the impurity effects on the SC 
on the surfaces of TIs. To construct the BdG equation, 
we perform the mean-field approximation 

C iTcr C jT' a' C jr' a' Cira 

{ ( HT<T < ^jr'(T l ) C 3'r'tT' c irtT + c Lct C ]t' ct' ( C ji"' 'a' ' c iTij) 



( C iT<T C jT' cr') (cjr'o-'CiTo-), 



(3.9) 



(3.10) 



which leads to the BCS mean-field Hamiltonian 

where * = (* e ^ h ) T = (c [c t ] T ) T is a Nambu spinor 
and 

Hmf -{ At -[H + H imp + ^i} T 



E A 



I 
/ 



(3.11) 



is the matrix clement. Here, fi is the Fermi energy, 
is a constant energy shift, and / is the identity ma- 
trix. The constant energy shift E& originates from the 
terms like (cc)(c^c^) in eq. (3.9) and calculated as E& = 

Eijaa'rr' 97J'(4ra c \r'a')( c 3T'a'C lT(T )- Since the contribu- 
tion of this term is just a constant energy shift, we are 
able to neglect this term in determining the mean-field 
self-consistently. In eq. (3.11), the pair potential is intro- 
duced as 



&Z>(h3) = 9ll'{i,j){c%TcrCjr'cT')- 



(3.12) 



The 4x4 pair potential matrix A is expressed as 



ITS 



EE A EH 



A 



EE 



.EH 



ti 

\HE A HH aHE \EH 

^tt ^tt °w ^ w 
AtV 3 AW A m B Af H 



it "4t 

HE AHH ^HE fcHH 



it 



,Eh 
Hi 

.Hi; 

hi 



(3.13) 



Note that A££, = — AJ/^.(j, «) is satisfied because of 
the Pauli exclusion principle. 

By diagonalizing Hmf, we obtain the excitation energy 
spectrum of the Bogoliubov quasiparticles E„ and the set 
of eigenvectors w v (it or) corresponding to amplitudes of 
the quasiparticle wave functions, that is, 



Hmfw v = E v w v . 



(3.14) 



The indices i, r, a, and r represent the site, subbands, 
spin, and particle-hole indices, respectively. 

The mean-field Hamiltonian has the particle-hole sym- 
metry (PHS), i.e. 



{C, Hmf} — 0, 



(3.15) 



where C = f x K and r x transforms a particle and a hole 
each other by the Pauli operator in the particle hole 
space and K is the complex conjugate operator. Then 
w v and Cw v are referred to as a particle-hole pair, i.e., 
HmfCw u = —EyCvju is satisfied. For convenience, we 
redefine the index v of an eigenvalue so that w v and W- v 
are a particle-hole pair, where v > and E v > 0. Then, 
w_ v = Cw v and E_ v = —E v hold. 

Then we introduce the creation and annihilation op- 
erators of quasiparticles, a„. Here, ^ and a v are related 
by the Bogoliubov transformation 



w v a v = (w u a v + w_ v a\) , (3.16) 



u>0 



and its inverse transformation 



w^, at 



(3.17) 



The relation a_„ = a J, holds because of PHS. Bogoliubov 
quasiparticles obey commutation relations of fermions so 
that {a v ,ajj} = 8 v>ll and {0^,0^} = are satisfied. The 
mean-field Hamiltonian is diagonalized as 



(3.18) 



u>0 



In order to calculate the pair potential, we intro- 
duce an imaginary time Green's function, G(r) = 
— (T r ^ r (r)^ ,T (0)). By the Bogoliubov transformation, the 
Green's function is expressed in the quasiparticle repre- 
sentation 



G(t) = ~y] w v wl (T t ol v (t)ol\) 



(3.19) 



The pair potential defined in eq. (3.12) is calculated 
from the anomalous part of the Green's function, i.e., 

= 5 -',(i,j)£{[l - f(E v )] Ul/ (iToK(jT'*') 



v>0 

fiE^vtiiTaJu^jr'a')}, 



(3.20) 
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where the two vectors u and v are introduced: 

u v {it<j) — w v {iT(je),v v {iTo) = w v (irah). (3.21) 

Here, u v and v v correspond to the amplitudes of the par- 
ticle and hole wave functions, respectively. In eq.(3.20), 
f(E) = (e l3E + is the Fermi-Dirac distribution func- 
tion. Note that u* = u_„ is satisfied because of PHS. 

3.2 Symmetry of the order parameter 

Wc first consider the impurity-free case and this sub- 
section is devoted to remarks satisfied in the absence of 
impurities. Wc analyze the order parameter, which is de- 
termined from the real-space BdG equation eq. (3.20). 
In this subsection, we concentrate on the order parame- 
ter A, which does not contain the interaction coefficents 
g^.,(i,j) introduced in eq. (3.7): 



result of the calculation is 



' v"> j) — (ciraCjT 'a' ) - 



(3.22) 



This order parameter A represents the supcrfluid density. 
The interaction parameter in eq. (3.7) is 

StTM = sifM = 9n H (i,i) = 9" H (i,i) = 2 eV; 

(3.23) 

for the other components, g r J a ,{i,i) = 0. By choosing the 
interaction coefficients as in eq. (3.23), we only consider 
onsitc attractive interactions between electrons in the 
same subband. We take the system size as N x = N y = 20 
and N z = 8 in this subsection. 

First, we show the onsite components A(i, i) of the 
order parameter obtained using the BdG equation. For 
the onsite component, the result of the calculation is 

/ o o m(z) im(z)\ 
o o -i m (z) m (z) 

-7?! (Z) i7 ?2 ( 2 ) 

\-m(z) -7/3(2) o o / 

1.45 x 10" 3 , 772(1) = 1-35 x 10~ 3 , 
1.05 x 10" 3 . (3.24) 



%(1) 



where 771 and 773 are even functions of z, while 772 is an 
odd function of z with the symmetry center at z = 4.5, 
i.e., the relations 771(2) = 771(9 — z), 772(2: — 9) = —7/2(9 — 
z), and 7/3(2: — 9) = 7/3(9 — 2) hold. The amplitude of 
each 77 decreases when 2 is nearer to the the symmetry 
center, for example, 771 (1) = 1.45 x 10~ 3 ,7/i(2) = 3.05 x 
10- 4 ,7/i(3) = 2.39 x 10~ 5 , 7/i(4) = 1.30 x 10~ 6 ; therefore 

7/ 1 (l)>77l(2)>7 7 l(3)>7/ 1 (4), 

holds. This amplitude dependences on 2 imply that the 
Cooper pairings are mainly formed from the helical Dirac 

electrons, which are localized at the surface. Since A(i, i) 
is an antisymmetric matrix, each onsite order parameter 

77 is a singlet component. Here, A(i,i) is independent of 
x and y because of translational symmetry. 

Then, we show a part of the off-site components 

A(i, i + e x ) of the order parameter obtained using the 

BdG equation. Here, A(i, i + e x ) represents the coherence 
between the two sites neighboring in the x-dircction. The 



f-Vi(z) iv' 2 ( z ) V' 3 (z) if/lM \ 

-Vaiz) ivi(z) -V'i( z ) -i^ 2 W 
\-iv' 4 (z) -77^(2) -irj 2 (z) rj' 5 (z) J 



7/i(l) = 1.63 x 10" 4 ,77 2 (1) = 1.71 x 10 



1.18 x 10" 3 , 7/4(1) 
1.60 x 10 _4 ,7/ 6 (1) 



1.02 x 10 -3 , 
1.12 x 10~ 3 , 



(3.25) 



where 7/2, rj' 3 , and r]' 6 are even functions of 2, while 77^ , 774, 
and 775 are odd functions of 2 with the symmetry center 
at 2 = 4.5. Since r)\_, r}' 2 , and rj' 5 are symmetric compo- 
nents of the above matrix, they are triplet components. 
In contrast, 773, 7/4, and ?7g are antisymmetric components 
of the matrix and are singlet components. 

The above results support the notion that an onsite 
s-wave interaction induces the order parameter with a 
mixture of singlet and triplet components, which agrees 
with the results of the idealistic helical Dirac electron 
model in §2. 

Generally, the order parameter is written as 

A(x, y, 2; x + t x e x ,y + t y e yi 2) 
= A(t, z)(t x I ®I-t y I® (icr 2 )) 
+ B{t, z)(t x T z ®I- t y r z ® (icr*)) 
+ C(t, z)(t x T x ® (-icr z ) - t y T x <X> 7) 
+ D(t, z)I ® (ia y ) + E(t, z)t z ® (ia y ) 



+ F{t,z)i 



(3.26) 



tl +ty. The above results are in the case of 



where t - 

t x = t y = and t x = l,t y = 0. Here, C, D, and E are 
even functions, while A, B, and F are odd functions of 2. 
In the above equation, a acts on the spin basis and the 
pseudospin Pauli matrix r acts on the basis of E and H 
subbands. 

According to the transformation rule in Appendix 1, 
the order parameter of the form in eq. (3.26) is invariant 
under a set of symmetry operations that belong to D4h- 
Since the noninteracting original Hamiltonian is invari- 
ant under the symmetry operations in D4h, the invari- 
ance of the order parameter under the operation indi- 
cates that the SC realized in the model breaks no addi- 
tional spatial symmetries. 

3.3 Stability against impurities 

Now, we analyze the impurity concentration and im- 
purity strength dependences of the pair potential for the 
tight-binding model of Bi2Sc3. Wc assume that impuri- 
ties are located on one of the two surfaces. We concen- 
trate on a pair potential on a surface with impurities. We 
define a surface with impurities as 2 = 1. 

By using the method of efficiently compensating for 
the change in the density of states (DOS) explained in 
§A.3 in the case of the reference systems, we analyze the 
impurity concentration dependence of the pair potential 
of Bi2Se3. Here, note that there are 16 components in 
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the BdG equation for Bi 2 Se 3 . Therefore, one might spec- 
ulate that the method of compensating for the change 
in the DOS for the BdG equation is complicated. How- 
ever, by choosing the Hamiltonian parameters used in 
§3.2, we can focus on and A^ . This is because we 
only consider onsite attractive interactions between elec- 
trons in the same subband. In the following discussion, 
we concentrate on the dependences of on the scat- 
tering strength and impurity concentration, and A^ is 
written as A. 

In order to describe the behavior of the pair potential 
purely due to relaxation processes, we concentrate on the 
value 

1 + 8A(u,nj) = ((A(u, ni )) x ) imp 2 ^ 
A (A (DOS) (u,n,i)) imp ' 

which represents the relative value of the pair potential 
due to the relaxation. In eq. (3.27), the quantities Ao, 
SA, A(u,rii), and A(dqs)( u ! n i) represent the impurity- 
free pair potential, the change in pair potential purely 
due to the relaxation, the result of the BdG calcula- 
tion, and the relaxation-ignored pair potential, respec- 
tively The definitions of these quantities are introduced 
in §A.3. Note that A(u, n{) depends on the site and impu- 
rity configuration, and A( DO S)( u 7 n i) depends on the im- 
purity configuration. Moreover, (• • • ) x and (• • • ); mp rep- 
resent the average over sites and impurity configuration, 
respectively, whose definitions are given in §A.3. 
In the model, we estimate A( D qs) as 

A(dos)(«, ni) = A 



A(DOS)(f,".) 



■sj E v (u,ni) 2 +A (DOS} (u,ni) : 



■■p(E u (u,rii),z = 1) 



An 



p<f>)(E 0v ,z=l) 



(3.28) 

where p(°'(E,z = 1) is the DOS at the surface with- 
out impurities, while p(E, z = 1) is that with impuri- 
ties. Here, Eq v and E„(u,rii) represent the energy spec- 
tra without and with impurities, respectively. We note 
that E v (u,n{) depends on the impurity configuration. 

We show the results of l + 5A(it, rn)/Ao in Figs. 3 and 
4. Figure 3 shows the impurity concentration dependence 
of the pair potential for the TRS impurities. In the figure, 
the results for the scattering strength of the choices u = 
0.1 and u = 20 are plotted. 

For u = 0.1, the inset in Fig. 3 shows that the reduc- 
tion in the pair potential is proportional to the square 
of the impurity concentration, i.e., cx nf for small n;. 
Since the order of the pair potential is 10 -3 , u = 0.1 is 
much larger than the pair potential. Therefore, it is be- 
yond the applicability of the AG theory. The robustness 
of the SC obtained in the BdG calculation indicates the 
stability of the SC beyond the pcrturbativc range. For 
u = 0.1, the pair potential vanishes at ~ 5%. On the 
other hand, when the strength of impurity potential is 
u = 20, which indicates a strong impurity potential re- 
garded practically as a lattice defect, the pair potential 
vanishes at rtj ~ 4%. 



In Fig. 4, we compare the results for TRS scattering 
with magnetic scattering. In the case of magnetic scat- 
tering, the pair potential decreases linearly for small con- 
centrations. We find that the reduction rate of the pair 
potential depends on the direction along x or z of the 
magnetization for magnetic impurities. We deduce that 
a difference between the two directions remains even at 
an impurity concentration range higher than 0.25% in 
spite of the large error bars. In fact, the difference in the 
impurity concentration dependences at 0.25% is reliable 
and statistically meaningful. Actually, according to the 
AG theory, the reduction rates of the pair potential are 
the same in both cases. Therefore, a difference in the 
impurity concentration dependence arises from higher- 
order perturbations. This difference is mainly caused by 
the difference in the change in DOS, which depends on 
the polarization. 33 



I<1 0.6 



+ 




Fig. 3. (Color online) Impurity concentration dependence of pair 
potential for TRS scattering. The abscissa represents the im- 
purity concentration rtj, while the ordinate represents 1 + 
SA(u, nj)/Arj. In the calculation, we fix the system size at 



N X N V 



400. The squares and circles illustrate the results for 



scatterings of strength u = 0.1 and u = 20, respectively. In the 
inset, the same quantity 1 + <5A(m, m)/Arj is plotted as a function 
of the square of the impurity concentration n? for u = 0.1. 



3-4 Spatial structure 

In this section, we analyze spatial structures induced 
by an impurity. We focus on the spatial structures of the 
pair potential and wave functions of Bogoliubov quasi- 
particles. To study the pair potential, we obtain a con- 
figuration of induced currents around an impurity. To 
study wave functions of quasiparticles, we analyze the 
properties of bound states around an impurity. Here, we 
assume that the impurity is located at (x, y) = (0, 0) and 
that the impurity potential u is strong. 

In the calculation, we choose the parameters of % as 



A, 



1 eV, A z = 0.5 eV, B x = 1.5 eV, B z = 0.3 eV, 



D X = D Z = 0, M = 0.5 eV, (3.29) 
<j££, (r, r) = 2 eV (onsite components), (3.30) 
g r J a ,{r, r') = 1 eV (r and r' are the nearest neighbors). 
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Fig. 4. (Color online) Impurity concentration dependence of pair 
potential for TRS and magnetic impurity scatterings. The ab- 
scissa represents the impurity concentration n;, while the ordi- 
nate represents 1 + SA(u, nj)/Ao. In the calculation, we fix the 
system size at N x N y = 400. The plus and cross illustrate the re- 
sults for the TRS impurity scattering when u = 0.1 and u = 20, 
respectively. The squares and circles illustrate the results for the 
magnetic impurity scattering, which were polarized along the al- 
and z-directions. 



(3.31) 

We take the system size as N x = N y = 20 and N z = 8. 

First, we show the spatial dependence of the onsite 
s-wave components around the impurity. Here, s-wave 
components arc defined as 

A£>) = A£'(r,r) 



A£ s (r) = AJf(r,r) 



(3.32) 



Figure 5 shows the spatial dependence of one of those 
components A^(r). The amplitude of A^(r) nearly 
vanishes at the impurity site. This is because the impu- 
rity potential is so strong that it is nearly regarded as a 
lattice defect and the electron density is extremely small 
at the impurity site. 



o 
i— i 

X 
< 




nian with an impurity potential. Because the system 
has the C4 rotational symmetry around the impurity 
site, the impurity site is a singular point of the vector 
field, (Re[A^(r)],Im[A^(r)]), where Ar^(r) is zero. 
Without the impurity, this component vanishes because 
the lattice translational symmetries along the x- and y- 
directions exist and every site is the rotational center. 

Figure 7 shows the spatial configurations of two 
normalized vector fields, F^. H (r)/\F^ H (r)\ and 
F^ H {r)/\F^ H {r)\, where 



U 



and 



■ tt 



(r) 



(Re[Af t >)],Im[Af t >)]), 



(3.33) 



F£» = (Re[Af*(r)],Im[A£>)]) . (3.34) 

Since we fix the length of the arrows representing the vec- 
tor fields, only the directions of the arrows are meaning- 
ful, which correspond to the complex phases of A^S(r) 
or A^p (r). At singular points, where \F\ = holds, the 
arrows are not plotted. In each component, there exist 
vortices at (0, 0), (0, 10), (10, 0), and (10, 10). In the A^ H 
component, the vortices at (0,0), (0, 10), and (10,0) are 
circulating clockwise while the vortex at (10, 10) is circu- 
lating counterclockwise. In contrast, the vortices of A^ H 
have opposite chiralities, i.e., 



^-jds-F^ H {r)/\F^ H {r)\ 



■-<t>ds-F*f(r)/\F*f(r)\ 



(3.35) 



In cq. (3.35), the chiralities of A^(r) and Af^f(r) 



F$ H (r)/\F{ 



2tt 



around a vortex are respectively defined as 

f t "(r)| and ± § ds- F™(r)/\F™(r 
an integral path around the vortex. 



f> ds ■ 
with 



Kit 

< 



Kit 
< 




1e-10 
1e-13 



Fig. 



5. Spatial dependence of |A^Fjp| 



Fie 



6. Spatial dependences of (a) |A^| and (b) In \A^\ 



Here, 



lA^f^l = holds at the origin and is not plotted in the lower 
panel. 



We then show the spatial dependence of the onsite p- 
wave components. Figure 6 shows A^S(r) = A^(r,r) 
which is one of the onsite p-wave components. This 
component reflects the existence of inter-orbital Cooper 
pairings. At the impurity site, this value is zero in 
agreement with the symmetry of the original Hamilto- 
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ds.F*f(r)/\F™(r)\ 



±-jds.F™{r)/\ F] 



EE 

44- 



(r)|, (3.40) 





x 



Fig. 7. (Color online) Spatial configuration of normalized vector 



fields, (a) 



and (b) Ffl H / 



F U, H \- The direction of 



the arrow represents the complex phase of A®^ 1 or A^H. At 
singular points, where |.F| = holds, the arrows are not plotted 



the induced circular electric current of each spin com- 
ponent flows in opposite directions. In eq.(3.40), the 
chiralities of Tj^_ and rif^L around a vortex are re- 
spectively defined as ^ § ds ■ FE^_{r)/\FR^_{r)\ and 
± § ds-F^i (r ) / 1 Ff^ (r)|, with an integral path around 
the vortex. These circulating currents represent spin cur- 
rents. 



In the present model, the off-site p-wave components 
A Px and A p exist simultaneously. They are calculated 
from the off-site pair potentials as 

)-X£(r,r-ei 



o 

, — i 

X 
63t 



(i = x,y). 




Then wc introduce the representation of p x ± ip y com- 
ponents as 

V Z'±(r) = \{K T ^{r) T iA^W), (3.37) 

where rj± corresponds to the p x ± ip y component. This 
representation is useful when describing the pair poten- 
tial of chiral p-wave-like SCs. 

Without impurities, (A^ p ^, A^ ) consists of 

components and (AJJ p ^, AJJ p ) consists of rj+ compo- 
nents. These results are consistent with the results in § 
3.1. For example, since A^ p (x,y,z = 1) = A(l,l) + 
B(l,l) and AfE )y (x,y,z= 1)"= -i[A(l, 1)+B(1, 1)] arc 
derived from eq. (3.26) without impurities, ij^^ vanishes 
and r)®®_ remains. 

We show the spatial structures of A^. E and A^P . Fig- 
ure 8 shows r/^_, which is the dominant component of 
A^R . This configuration is the same as that of r]^_. 
Figure 9 shows the absolute values of induced by 
the impurity. From the spatial configuration on the log- 
arithmic scale, we find that there are six vortices in the 
system. Figure 10 shows the spatial configurations of two 
normalized vector fields, FRf/\FRf\ and Ff^/\F^_\, 
where 



(3.36) 

Fig. 8. Spatial dependence of [fj^rj, which is the dominant com- 
ponent of A®® . This configuration is the same as | r) | . 
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and 



(r) = (Re[ryf4(r)],Im[r,f4(r)]) , (3.38) 



(3.39) 



Fj^(r) = (Re[^(r)],Imfo££(r)]) 

Since we fix the length of the arrows representing the vec- 
tor fields, only the directions of the arrows are meaning- 
ful, which correspond to the complex phases of (r) 
or Tffu- ( r )- At singular points, the arrows are not plotted 
as in Fig. 7. Since 17^. and flf^L have opposite chiralities, 



Fig. 9. Spatial dependences of (a) and (b) In | t?^5[^_ | - These 

configurations are the same as and In respectively. 



In this subsection, we observe the spatial structures 
of typical pair-potential components. For the complex 
phase of s-wave and p-wave components, qualitatively 
different spatial configurations arc obtained, both of 
which are allowed in the original symmetry of the system. 
In particular, the complex phase of p-wave components 
around impurities should be observed as spin currents 
around them. 
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10. (Color online) Spatial configurations of normalized vec- 
tor fields, (a) F$E/\ Fj EE\ and ( b ) F^/\F^_\. The direction 

of the arrow represents the complex phase of TO-S. or vf±L ■ At 
singular points, where |.F| = holds, the arrows are not plotted. 
The two opposite components have opposite chiralities. 



3. 5 Wave function of Bogoliubov quasi-particle 

Now, we study the energy spectra and wave functions 
of a Bogoliubov quasiparticle composed of SSs. In the 
present calculation, there are eight states near the Fermi 
energy Without impurities these eight states are degen- 
erate since the system has a four-fold rotational symme- 
try and an inversion symmetry. By introducing an im- 
purity potential, these states split. In Fig. 11, we show 
the ratio of the energy level splitting E imp to the BCS 
energy gap Ebcs caused by an impurity when the sys- 
tem size changes. Since i?i mp / -Ebcs is scaled to zero in 
the thermodynamic limit, it supports the notion that the 
impurity level does not split off from the gap edge. This 
means that the state docs not develop into a mid-gap 
bound state. 




0.02 0.04 0.06 0.08 0.1 0.12 0.14 

Fig. 11. Relationship between system size and ratio of energy 
level splitting E- lulp to BCS energy gap -Ebcs caused by an im- 
purity. Inset shows a log-log plot for the same data. 



In Fig. 12, the amplitude of the quasiparticle wave 
function that has an energy nearest to the Fermi energy 
is shown. This wave function is spatially extended and 
supports the notion that, in the thermodynamic limit, 
the state in Fig. 11 docs not seem to reduce to a mid-gap 
bound state, in agreement with the absence of a bound 
state. 




10 



Fig. 12. Amplitude of the quasiparticle wave function for the 
state whose energy is the nearest to the Fermi energy. 



If there were an energy level of a mid-gap or gapless 
state in a SC, it would not be able to penetrate the bulk 
due to the existence of a SC gap. Therefore, such a mid- 
gap state has to be bounded near the surfaces or edges 
of the SC. The mid-gap state is referred to as an An- 
dreev bound state (ABS). Because the SC is a full gap 
SC in the 2D surface, the stability of ABSs is related to 
the topology of the superconducting gap; if the super- 
conducting gap is topologically nontrivial, ABSs should 
be topologically protected and stable, and such a SC is 
called a topological SC, which is analogous in stability 
of SSs in TIs. 

Our result supports the notion that there are no ABSs 
that are mid-gap states formed around the impurity. The 
absence of ABSs represents the notion that the supercon- 
ductivity gap induced by an s-wave attractive interaction 
on the surfaces of TIs is topologically trivial, i.e., the 
present SC is not a topological SC. On the other hand, 
spinless chiral p-wave SCs, which have similar forms of 
pair potential to the present SCs, are topological SCs. 14 
This means that SCs on the surfaces of TIs and spinless 
chiral p-wave SCs are topologically different. 

4. Conclusion and Discussion 

In §2, we have shown that unconventional SCs induced 
by the s-wave attractive interaction on the surfaces of 
TIs are robust against TRS disorders in idealistic Dirac 
electron models. This is because of the cancellation of 
two phase factors, i.e., one from the pairing potential and 
the other from the scattering factor of Dirac electrons. 
In contrast, unconventional SCs studied in the literature, 
such as d-w&ve and chiral p-wave SCs, 23 are sensitively 
suppressed through scattering by a tiny concentration 
of impurities because of the phase factor of the pairing 
potential. 

In §3, by numerically analyzing the Bogoliubov de- 
Gennes equation beyond the perturbative regime for im- 
purities, we have obtained the result that SCs induced by 
the s-wave attractive interaction on the surfaces of TIs 
are stable against TRS impurities since the pair potential 
does not decrease linearly within a range of small con- 
centration. This result is consistent with that discussed 
in §2. Moreover, we have found that the robustness is 
observed even when the impurity potential is strong be- 
yond the perturbation theory. We have also found that 
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the reduction rate depends on the spin polarization of 
magnetic impurities. This is because the change in the 
DOS depends on the polarization. 33 This implies that 
the superconducting gap in the present SC is topologi- 
cally trivial, i.e., the present SC is a trivial SC. Generally, 
the existence of ABSs and the stability of SCs are closely 
related, i.e., SCs that are stable against TRS impurities 
have ABSs around impurities. According to our results, 
this relation appears to be satisfied in the present SC on 
the surfaces of TIs. 

Finally, we describe issues left for future study. In this 
paper, we have studied the superconductivity stabilized 
by an internal s-wave interaction in order to clarify fun- 
damental impurity effects. Recently, however, supercon- 
ductivity introduced by the proximity effect has also ac- 
tively been studied. 12 ' 17 Impurity effects on the SC in- 
troduced by the proximity effect are also intriguing. 

Moreover, stability against vortices is important. This 
is because, in realizing a Majorana bound state around 
a vortex, 12 the backscattering of Dirac electrons due to 
a magnetic field, which breaks TRS, can be destructive 
in terms of the stability of SSs. 
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Appendix: 

A.l Symmetry operation 

Here, we discuss symmetry operations that preserve 
spatial symmetries that the Hamiltonian in eq. (3.2) has. 
The symmetries of superconducting pair potentials are 
determined by irreducible representations of a group con- 
sisting of these symmetry operators. 

First, we introduce a matrix representation of sym- 
metry operators as follows: A symmetry operator on a 
Hilbert space is denoted as V, while its matrix represen- 
tation acting on creation (annihilation) operator-vectors 
P (P^) is defined as 

Vc^V- 1 = JP, (A-l) 

PcP- 1 = P^c, (A-2) 

where P depends on the choice of the basis. 

Then the BCS mean-field Hamiltonian is transformed 
under a symmetry operation P as 



pup- 




il PUP' 1 = H , i.e., PH P f = H and PAP T = A, 
is satisfied, the noninteracting Hamiltonian H$ and the 
order parameter matrix A are both invariant under the 
operation P. 



The effective tight-binding Hamiltonian Ho of Bi2Se3 
in §3 is invariant under the symmetry operations belong- 
ing to the point group D4h- These symmetry operations 
are characterized by the character table given in Table 
A.l. 

With the choice of the basis employed in this paper, 
we obtain matrix representation for each symmetry op- 
eration in D4h, as shown in Tables A. 2. and A. 3. Table 
A. 2 corresponds to the matrix acting on the spin indices 
and subbands indices E and H , and Table A. 3 corre- 
sponds to the matrix acting on spatial coordinates. In 
Table A. 2, we note that a acts on spin indices and f 
acts on subbands indices. When we assume an s-wave 
attractive interaction in the effective model, we obtain 
pair potentials belonging to the irreducible representa- 
tion Aij,.. 
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Tabic A-l. Character table for D 4 h point group. 





E 




c? 


2Co 


2(%' 


J 




w fl 


2a v 


2(Tw 
^ Jiy a 


Example of basis 


Ai„ 

rg 


1 


1 




1 


1 


1 


1 


l 


1 


1 


z 2 or x 2 + y 2 


^g 


1 


1 


J 


-1 


-1 


1 


1 


l 


-1 


-1 


xv(x 2 — v 2 ) 


Big 


1 


-1 




1 


-1 


1 


-1 


1 


1 


-1 


x 2 — ?/ 2 


1 V)., 


1 


-1 




-1 


1 


1 


-1 


l 


-1 


1 


a/ y 


Eg 


2 





—2 





o 


2 





-2 





o 


\ — ZX Z1l\ 


Ai„ 


1 


1 




1 


1 


-1 


-1 


-1 


-1 


-1 


xaiz( X 2 — 1l 2 \ 


A 2u 


1 


1 




-1 


-1 


-1 


-1 


-1 


1 


1 


Z 




1 


-1 




1 


-1 


-1 


1 


-1 


-1 


1 


xyz 


B2u 


1 


-1 




-1 


1 


-1 


1 


-1 


1 


-1 


z(x 2 - y 2 ) 


E u 


2 





-2 








-2 





2 








{x,y} 
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Table A-2. Table of transformation matrices for spin basis and 
subbands indices. 





(E,H)®(t,l) 


E 




I® I 




2C 4 


75 / 


g)(/ + icr z ), 




c\ 


il ® cr z 




2C' 2 




U <g> <7 X , U ig) (Ty 




2C» 


75 /g 


) ((Tit + CTj,), -^1 ® (o-j, - 


- °v) 


I 


T z ® / 




2/C4 = 2S 4 


75- 


8 (7 + icr z ), -l=r z ® (/- 


- ic z ) 


°h 




ir z <g> <r z 




2(T V 




ir z ® (Ta; , ir z ® <r H 




2<r d 


75 T ^ 


) (fas + -ijr z g) ((Ta, 





Table A-3. Transformation for spatial coordinates. 





(x, y, 2) 


E 


{x,y,z) 




2C 4 


(y, ~x,z), (-y,x 


z) 




(-x, -y,z) 




2C^ 


(x, -y, -z), {-x,y, 


-z) 


2C» 


(y,x,-z), (-y,-x 


-z) 


I 


(— a?> -J/! -2) 




2IC4, = 2S 4 


(-!/,*, -2), (2/; -x 


-') 


Oh 


(x,y, -z) 




2a v 


(-x,y,z), (x,-y, 


z) 


2a d 


(-y, -x,z), (y,x 


*) 



A. 2 Reference systems 

When we analyze impurity effects on 2D superconduc- 
tivities specific to those arising in the surface states of 
3D TIs, we need to study those on the 2D SC in a topo- 
logically trivial reference system. As a simple reference 
system, here we introduce a 2D tight-binding model on 
a square lattice. 

The Hamiltonian of the 2D tight-binding model with 
an s-wave attractive interaction is defined as 

ft = 1 4x c ^-sXXW c !i c 4 

Ni 

i— 1 ia 

where t, g, u, and fj. are the hopping between nearest- 
neighbor sites, the amplitude of the on-site attractive 
interaction, the impurity potential, and the chemical po- 
tential of the system, respectively. Here, n(i) is the loca- 
tion of the i-th impurity and JVj is the number of impu- 
rities. When we consider magnetic scattering (polarized 
along the z-axis), we substitute 

AT, 

XX(i)t C ™Wt " C i(zW C ««l)> ( A ' 5 ) 
i=l 

for the impurity term in eq.(A-4), 

U H C l W a C «W-- ( A ' 6 ) 
t=l 



Author Name 15 

Similarly to the effective model for Bi2Se3, we intro- 
duce a mean-field decoupling with real-space order pa- 
rameters Aj being self-consistently defined by 

A. t = g(c it cn). (A-7) 

Note that, when u — 0, the order parameter Aj is homo- 
geneous and obtained using 

A, = A , 

A. 3 Method of compensating for change in density of 
states 

In this subsection, we introduce a simple scheme for 
reducing finite-size effects in impurity effects on SC. A 
dominant finite size effect comes from the change in the 
DOS due to the introduced impurity potentials, which 
is assumed to be negligible in a thermodynamic limit. 
Then, we test the validity of the method by analyzing 
numerical results on the s-wave SC in the tight-binding 
model on the square lattice, introduced in the previous 
subsection, in comparison with the AG theory 

The change in the DOS is assigned as a higher-order ef- 
fect and neglected in the AG theory for continuum mod- 
els. However, in finite-size BdG calculations, the change 
in the DOS arising from the impurities quantitatively af- 
fects the order parameter, in addition to relaxation or 
pair-breaking processes due to impurities, while this is 
assumed to be negligible. Therefore, we have to com- 
pensate for this change when one wishes to estimate the 
genuine reduction in the order parameter purely arising 
from impurity relaxation processes. Here we remind the 
readers that the change in the DOS and the effects of the 
impurity relaxation processes correspond to the diagram 
in Figs. 1(a) and 1(b), respectively. 

When we neglect such a change in the DOS, the pair 
potential with impurities is calculated in the perturba- 
tion regime as 

A (It, Hi) = A + SA{u,m), (A-9) 

6A ^ = -^ry (a-io) 

according to the AG theory 24, 30 where u and ri; is the 
strength of impurity potentials and the impurity con- 
centration, respectively. In eq.(A-9), Ao is the impurity- 
free pair potential. In eq.(A-lO), t s (u, n\) is the impurity 
relaxation time contributing to the reduction in the or- 
der parameter. This relaxation time is determined by the 
strength of impurity potentials u and the impurity con- 
centration rij as well as the symmetries of the impurity 
scattering and order parameters. 23,24 

However, in numerical solutions of the BdG equation 
on finite-size systems, the reduction in the order param- 
eter is not fully given by SA(u, nj) in eq.(A-9). We need 
to take into account the changes in the DOS as well. In 
order to subtract the reduction due to the changes in 
the DOS, we introduce a "relaxation-ignored" pair po- 
tential Amos); m which only the change in the DOS by 
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the impurities is taken into account, while the impurity 
relaxation times are neglected. By using the relaxation- 
ignored pair potential A( DO S)j eq.(A-9) is replaced by 

A(tt,7ii) = A (dos) (m, m) + 5A(u,ni), (A-ll) 



when 



u 2 riiN 
A(u,rn) 



< 1 



(A-12) 



holds, where No is the DOS at the Fermi energy. Here, 
u 2 niiVo, the numerator on the left-hand side of eq.(A-12), 
is the same order as SA{u,rii). 

Then, we introduce an equation by which the 
relaxation-ignored pair potential is calculated and ver- 
ify the validity of the method by analyzing s-wave SC 
in a tight-binding model for the square lattice. In the 
regime of small impurity concentration, we define the 
relaxation-ignored pair potential Apos) horn the self- 
consistent equation 



A(DOS)0,ni) = f 



^(DOS) 



(u,m) 



S A^/E v (u,mY + A ( Dos)(w,«i) 

(A-13) 

where E v (u, n{) is the set of eigenvalues when the sys- 
tem has no attractive interaction but has impurities. By 
the estimation using eq.(A-13), we consider the change in 
the order parameter due to the shift in the energy spec- 
trum. The relaxation processes are not contained in the 
estimate of A(dos)> because they appear as shifts in the 
imaginary parts of the quasiparticle energies. 

The estimation of the relaxation- ignored pair potential 
by eq.(A-13) enables us to compensate for the change in 
the DOS in the BdG calculations. When we focus on the 
relative reduction in the pair potential purely from the 
impurity relaxation processes, we should concentrate on 
the quantity 5A(u,rii)/Ao, which represents the reduc- 
tion. By using eq.(A-ll), we introduce an equation that 
associates the results of the BdG calculations with the 
reduction of the pair potential purely due to the relax- 
ation processes. The equation is 

5A(u,rii) _ ((A(u,ni)) x ) imp 



1 



(A(DOS)(«,«-i))i 



(A-14) 



which is valid in the ranges of small n\ and u. Since the 
pair potentials obtained by the BdG calculations have 
spatial and impurity-configuration dependences, we take 
two types of averages: (• • • ) x and (• • • )i mp . The average 
(■ ■ - ) x means the spatial average, i.e., 



1 S 



(A-15) 



where S is the system size and Ai is a quantity depending 
on the site i. The other average (• • • ); mp represents the 
average over impurity configurations, i.e., 



Ay 



where N c is the number of impurity configurations and 



Bj is a quantity depending on the impurity configuration 
j. Moreover, note that we have to take Ados after the 
impurity-configuration average, since it also depends on 
the impurity configuration. 

Then we can compare the results of the BdG with 
those obtained using the AG theory. By using the AG 
theory, the left side of cq.(A-14) is estimated as 



1 



<5A(it, rii) 



= f - 



(A-17) 



A 4t s (m, n ; )A 

In the cases of the TRS and magnetic scatterings, the 
relaxation times r s are obtained using 

1 



(TRS) 



and 



(mag) , 



l-KTiiU N n , 



(A-18) 



(A-19) 



respectively, where Nq is the DOS at the Fermi energy 
without impurities, i.e., 



2irS 



(A-20) 



Here, note that, in the finite size system, the DOS is not 
well-defined. However, we can introduce a reasonable es- 
timation of the DOS Nq by reconsidering how the DOS 
Nq appears in the AG theory. In the AG theory, the quan- 
tity Nq is replaced by the the DOS in the thermodynamic 
limit 



1 



2ttS 



dw 



a.' 



uj ui + e 



(A-21) 



lie 



where w is the energy cutoff. We adopt uj as the band- 
width, i.e., ui = 8, and estimate the DOS as N = N^ 
even in finite-size systems. 

By comparing the two estimate of 1 + <5A ^"'- > by eqs. 
(A-13) and (A-14), we show the validity of our method of 
compensating for the change in the DOS. Figures A-l and 
A-2 show the scattering strength and impurity concentra- 
tion dependences of the pair potential, respectively. By 
the analyses of both dependences, we find an agreement 
between the two different approaches: the AG theory and 
real-space BdG calculation. Therefore, we conclude that 
our method of compensating for the change in the DOS 
is valid for analyzing the impurity concentration depen- 
dence of the order parameter. 
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Fig. A-l. (Color online) Scattering strength dependence of pair 
potential for magnetic scattering. The abscissa represents the 
square of the strength of the impurity potential u, while the 
ordinate represents 1 + 8A(u, ri;)/Ao, which corresponds to a 
relative change in the pair potential. We take the system size as 



N x 



N v 



solid line is a theoretical estimate by the AG theory. The circle 
symbol corresponds to the solution of the BdG equations. 
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Fig. A-2. (Color online) Impurity concentration dependence of 
the pair potential for TRS and magnetic scatterings. The ab- 
scissa represents the impurity concentration n;, while the ordi- 
nate represents 1 + 8A(u, ni)/Ao. In the calculation for the TRS 
scattering, we fix the scattering strength to u = 0.1 and the num- 
ber of impurities to Ni = 1. We alter the impurity concentration 
by changing the system size N x N y . In the calculation for the 
magnetic scattering, we fix the scattering strength tou = 0.1 and 
the system size to N x N y = 400. The solid and dashed lines are 
the theoretical estimates for the TRS and magnetic scatterings 
by the AG theory, respectively. The circle and triangle symbols 
correspond to the results of the TRS and magnetic scatterings 
by the BdG calculations, respectively. 
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